# Respectability politics and straight support for LGB rights
# Replication code
# Phil Jones
# 07.10.21

library(car)
library(survey)
library(ggplot2)
library(cowplot)
library(patchwork)
library(scales)
library(srvyr)
library(texreg)
library(poliscidata)
library(FindIt)


load("RespectabilityPolitics_Study1.RData")
load("RespectabilityPolitics_Study2.RData")

# Study1.s: weighted CCES data, straight respondents only
# Study2.s: unweighted Qualtrics data, straight respondents only 

phil.ci <- function(x){
  cis <- confint(x)
  ans <- cbind(as.matrix(x), cis)
  ans <- round(ans, digits=3)
  ans <- as.data.frame(ans)
  names(ans) <- c("est", "lo", "hi")
  ans
}



##########################################################################################
# STUDY 1: ATEs OF BEING IN A TWO-PERSON RELATIONSHIP
##########################################################################################

# Regression models
reg_case <- svyglm(case.outcome~treat, Study1.s)
reg_siml <- svyglm(similar.gaymen~treat, Study1.s)
reg_angr <- svyglm(felt.angry~treat, Study1.s)
reg_disg <- svyglm(felt.disgust~treat, Study1.s)
reg_prod <- svyglm(felt.proud~treat, Study1.s)
reg_serv <- svyglm(lgb.services~treat, Study1.s)
reg_rght <- svyglm(lgb.rights~treat, Study1.s)

# Estimate ATEs
d_mono <- data.frame(treat=0) 
d_poly <- data.frame(treat=1) 
fd_case <- phil.ci(predict(reg_case, d_poly, type="response")-predict(reg_case, d_mono, type="response"))
fd_siml <- phil.ci(predict(reg_siml, d_poly, type="response")-predict(reg_siml, d_mono, type="response"))
fd_angr <- phil.ci(predict(reg_angr, d_poly, type="response")-predict(reg_angr, d_mono, type="response"))
fd_disg <- phil.ci(predict(reg_disg, d_poly, type="response")-predict(reg_disg, d_mono, type="response"))
fd_prod <- phil.ci(predict(reg_prod, d_poly, type="response")-predict(reg_prod, d_mono, type="response"))
fd_serv <- phil.ci(predict(reg_serv, d_poly, type="response")-predict(reg_serv, d_mono, type="response"))
fd_rght <- phil.ci(predict(reg_rght, d_poly, type="response")-predict(reg_rght, d_mono, type="response"))
fd_study1 <- rbind(fd_siml, fd_case, fd_serv, fd_rght, fd_angr, fd_disg, fd_prod)

# Plot ATEs
row.names(fd_study1) <- 1:7
fd_study1 <- as.data.frame(fd_study1)
fd_study1$dv <- c("Felt similar to the gay men", "Support the gay men's case", "Support required service","Support LGB rights", "Felt angry", "Felt disgusted", "Felt proud")
fd_study1$dv <- factor(fd_study1$dv)
fd_study1$dv <- factor(fd_study1$dv, levels(fd_study1$dv)[c(3,2,1,5,6,7,4)]) 
fd_study1$grp <- c(rep("p1",2), rep("p2", 2), rep("p3",3))

Figure1 <- ggplot(fd_study1, aes(x=dv, y=est, ymin=lo, ymax=hi)) +
  geom_hline(yintercept=0, color="gray55", linetype=1, size=.5) +
  geom_errorbar(width=.1, size=.55)+
  geom_point(size=3, shape=21, fill="gray30") +
  scale_y_continuous(limits=c(-.3,.3), breaks=c(-.3, -.2, -.1, 0,.1,.2,.3), minor_breaks=NULL) + 
  ylab("\nAverage Treatment Effect of showing gay men in two-person relationship") + xlab("") +
  coord_flip() +
  theme(axis.title=element_text(size=11),
        axis.text.x=element_text(size=9, color="black"),
        axis.text.y=element_text(size=11, color="black"),
        axis.ticks=element_line(size=.5, color="gray55"),
        legend.text=element_text(size=11),
        legend.position="top")
Figure1 + facet_grid(grp~., scales="free", space="free_y") + theme(strip.text=element_blank())
ggsave("Figure1.pdf", width=7.5, height=6.5)




##########################################################################################
# STUDY 1: Do respondents’ predispositions moderate these effects?
##########################################################################################

# Models to estimate heterogeneous treatment effects on three DVs

# LGB THERMOMETER
reg_case_therm <- svyglm(case.outcome~treat*lgb.therm, Study1.s)
reg_sim_therm <- svyglm(similar.gaymen~treat*lgb.therm, Study1.s)
reg_serv_therm <- svyglm(lgb.services~treat*lgb.therm, Study1.s)
# IDEOLOGY
reg_case_ideo <- svyglm(case.outcome~treat*ideo, Study1.s)
reg_sim_ideo <- svyglm(similar.gaymen~treat*ideo, Study1.s)
reg_serv_ideo <- svyglm(lgb.services~treat*ideo, Study1.s)
# PID7
reg_case_pid7 <- svyglm(case.outcome~treat*pid7, Study1.s)
reg_sim_pid7 <- svyglm(similar.gaymen~treat*pid7, Study1.s)
reg_serv_pid7 <- svyglm(lgb.services~treat*pid7, Study1.s)

# Estimate effects by each moderator
d_mono <- data.frame(treat=0, lgb.therm=seq(from=0,to=1,by=.1), ideo=seq(from=0,to=1,by=.1), pid7=seq(from=0,to=1,by=.1)) 
d_poly <- data.frame(treat=1, lgb.therm=seq(from=0,to=1,by=.1), ideo=seq(from=0,to=1,by=.1), pid7=seq(from=0,to=1,by=.1)) 
fd_case_therm <- phil.ci(predict(reg_case_therm, d_poly, type="response")-predict(reg_case_therm, d_mono, type="response"))
fd_case_therm$dv <- "case.outcome"
fd_case_therm$x <- seq(from=0,to=1,by=.1)
fd_case_therm$iv <- "lgb.therm"
fd_sim_therm <- phil.ci(predict(reg_sim_therm, d_poly, type="response")-predict(reg_sim_therm, d_mono, type="response"))
fd_sim_therm$dv <- "similar.gaymen"
fd_sim_therm$x <- seq(from=0,to=1,by=.1)
fd_sim_therm$iv <- "lgb.therm"
fd_serv_therm <- phil.ci(predict(reg_serv_therm, d_poly, type="response")-predict(reg_serv_therm, d_mono, type="response"))
fd_serv_therm$dv <- "lgb.services"
fd_serv_therm$x <- seq(from=0,to=1,by=.1)
fd_serv_therm$iv <- "lgb.therm"

fd_case_ideo <- phil.ci(predict(reg_case_ideo, d_poly, type="response")-predict(reg_case_ideo, d_mono, type="response"))
fd_case_ideo$dv <- "case.outcome"
fd_case_ideo$x <- seq(from=0,to=1,by=.1)
fd_case_ideo$iv <- "ideo"
fd_sim_ideo <- phil.ci(predict(reg_sim_ideo, d_poly, type="response")-predict(reg_sim_ideo, d_mono, type="response"))
fd_sim_ideo$dv <- "similar.gaymen"
fd_sim_ideo$x <- seq(from=0,to=1,by=.1)
fd_sim_ideo$iv <- "ideo"
fd_serv_ideo <- phil.ci(predict(reg_serv_ideo, d_poly, type="response")-predict(reg_serv_ideo, d_mono, type="response"))
fd_serv_ideo$dv <- "lgb.services"
fd_serv_ideo$x <- seq(from=0,to=1,by=.1)
fd_serv_ideo$iv <- "ideo"

fd_case_pid7 <- phil.ci(predict(reg_case_pid7, d_poly, type="response")-predict(reg_case_pid7, d_mono, type="response"))
fd_case_pid7$dv <- "case.outcome"
fd_case_pid7$x <- seq(from=0,to=1,by=.1)
fd_case_pid7$iv <- "pid7"
fd_sim_pid7 <- phil.ci(predict(reg_sim_pid7, d_poly, type="response")-predict(reg_sim_pid7, d_mono, type="response"))
fd_sim_pid7$dv <- "similar.gaymen"
fd_sim_pid7$x <- seq(from=0,to=1,by=.1)
fd_sim_pid7$iv <- "pid7"
fd_serv_pid7 <- phil.ci(predict(reg_serv_pid7, d_poly, type="response")-predict(reg_serv_pid7, d_mono, type="response"))
fd_serv_pid7$dv <- "lgb.services"
fd_serv_pid7$x <- seq(from=0,to=1,by=.1)
fd_serv_pid7$iv <- "pid7"


# Plot effects

# DV: Feeling similar to gay man
Figure2a <- ggplot(fd_sim_therm, aes(x=x, y=est, ymin=lo, ymax=hi, group=iv)) +
  geom_hline(yintercept=0, color="gray55", linetype=1, size=.5) +
  geom_ribbon(alpha=.2) +
  geom_line(size=1) +
  scale_y_continuous(limits=c(-.3,.3), breaks=c(-.3, -.2, -.1, 0,.1,.2,.3), minor_breaks=NULL) +
  scale_x_continuous(limits=c(0,1), breaks=c(0,1), labels=c("Most\ncold","Most\nwarm"), minor_breaks=c(.25,.5,.75)) +
  ylab("") + xlab("") + ggtitle("by LG thermometer rating") +
  theme(axis.text.y=element_text(color="black", size=9),
        axis.text.x=element_text(color="black", size=10, hjust=c(0,1)),
        plot.title=element_text(color="black", size=13, hjust=.5),
        axis.ticks=element_line(size=.5, color="gray55"))
Figure2a <- Figure2a + coord_cartesian(clip = "off") + draw_label("Effect of two-person relationship", fontface='plain', size=13, angle=90, x=-.23, y=0)
Figure2b <- ggplot(fd_sim_ideo, aes(x=x, y=est, ymin=lo, ymax=hi, group=iv)) +
  geom_hline(yintercept=0, color="gray55", linetype=1, size=.5) +
  geom_ribbon(alpha=.2) +
  geom_line(size=1) +
  scale_y_continuous(limits=c(-.3,.3), breaks=c(-.3, -.2, -.1, 0,.1,.2,.3), minor_breaks=NULL) +
  scale_x_continuous(limits=c(0,1), breaks=c(0,1), labels=c("Very\nliberal","Very\nconservative"), minor_breaks=c(.25,.5,.75)) +
  ylab("") + xlab("") + ggtitle("by ideology") +
  theme(axis.text.y=element_text(color="black", size=9), axis.ticks.y=element_blank(),
        axis.text.x=element_text(color="black", size=10, hjust=c(0,1)),
        plot.title=element_text(color="black", size=13, hjust=.5),
        axis.ticks=element_line(size=.5, color="gray55"))
Figure2c <- ggplot(fd_sim_pid7, aes(x=x, y=est, ymin=lo, ymax=hi, group=iv)) +
  geom_hline(yintercept=0, color="gray55", linetype=1, size=.5) +
  geom_ribbon(alpha=.2) +
  geom_line(size=1) +
  scale_y_continuous(limits=c(-.3,.3), breaks=c(-.3, -.2, -.1, 0,.1,.2,.3), minor_breaks=NULL) +
  scale_x_continuous(limits=c(0,1), breaks=c(0,1), labels=c("Strong\nDemocrat","Strong\nRepublican"), minor_breaks=c(.25,.5,.75)) +
  ylab("") + xlab("") + ggtitle("by party identity") +
  theme(axis.text.y=element_text(color="black", size=9), axis.ticks.y=element_blank(),
        axis.text.x=element_text(color="black", size=10, hjust=c(0,1)),
        plot.title=element_text(color="black", size=13, hjust=.5),
        axis.ticks=element_line(size=.5, color="gray55"))


# DV: Supporting gay man in case
Figure2d <- ggplot(fd_case_therm, aes(x=x, y=est, ymin=lo, ymax=hi, group=iv)) +
  geom_hline(yintercept=0, color="gray55", linetype=1, size=.5) +
  geom_ribbon(alpha=.2) +
  geom_line(size=1) +
  scale_y_continuous(limits=c(-.3,.3), breaks=c(-.3, -.2, -.1, 0,.1,.2,.3), minor_breaks=NULL) +
  scale_x_continuous(limits=c(0,1), breaks=c(0,1), labels=c("Most\ncold","Most\nwarm"), minor_breaks=c(.25,.5,.75)) +
  ylab("") + xlab("") + ggtitle("by LG thermometer rating") +
  theme(axis.text.y=element_text(color="black", size=9),
        axis.text.x=element_text(color="black", size=10, hjust=c(0,1)),
        plot.title=element_text(color="black", size=13, hjust=.5),
        axis.ticks=element_line(size=.5, color="gray55"))
Figure2d <- Figure2d + coord_cartesian(clip = "off") + draw_label("Effect of two-person relationship", fontface='plain', size=13, angle=90, x=-.23, y=0)
Figure2e <- ggplot(fd_case_ideo, aes(x=x, y=est, ymin=lo, ymax=hi, group=iv)) +
  geom_hline(yintercept=0, color="gray55", linetype=1, size=.5) +
  geom_ribbon(alpha=.2) +
  geom_line(size=1) +
  scale_y_continuous(limits=c(-.3,.3), breaks=c(-.3, -.2, -.1, 0,.1,.2,.3), minor_breaks=NULL) +
  scale_x_continuous(limits=c(0,1), breaks=c(0,1), labels=c("Very\nliberal","Very\nconservative"), minor_breaks=c(.25,.5,.75)) +
  ylab("") + xlab("") + ggtitle("by ideology") +
  theme(axis.text.y=element_text(color="black", size=9), axis.ticks.y=element_blank(),
        axis.text.x=element_text(color="black", size=10, hjust=c(0,1)),
        plot.title=element_text(color="black", size=13, hjust=.5),
        axis.ticks=element_line(size=.5, color="gray55"))
Figure2f <- ggplot(fd_case_pid7, aes(x=x, y=est, ymin=lo, ymax=hi, group=iv)) +
  geom_hline(yintercept=0, color="gray55", linetype=1, size=.5) +
  geom_ribbon(alpha=.2) +
  geom_line(size=1) +
  scale_y_continuous(limits=c(-.3,.3), breaks=c(-.3, -.2, -.1, 0,.1,.2,.3), minor_breaks=NULL) +
  scale_x_continuous(limits=c(0,1), breaks=c(0,1), labels=c("Strong\nDemocrat","Strong\nRepublican"), minor_breaks=c(.25,.5,.75)) +
  ylab("") + xlab("") + ggtitle("by party identity") +
  theme(axis.text.y=element_text(color="black", size=9), axis.ticks.y=element_blank(),
        axis.text.x=element_text(color="black", size=10, hjust=c(0,1)),
        plot.title=element_text(color="black", size=13, hjust=.5),
        axis.ticks=element_line(size=.5, color="gray55"))


# DV: Supporting required service
Figure2g <- ggplot(fd_serv_therm, aes(x=x, y=est, ymin=lo, ymax=hi, group=iv)) +
  geom_hline(yintercept=0, color="gray55", linetype=1, size=.5) +
  geom_ribbon(alpha=.2) +
  geom_line(size=1) +
  scale_y_continuous(limits=c(-.3,.3), breaks=c(-.3, -.2, -.1, 0,.1,.2,.3), minor_breaks=NULL) +
  scale_x_continuous(limits=c(0,1), breaks=c(0,1), labels=c("Most\ncold","Most\nwarm"), minor_breaks=c(.25,.5,.75)) +
  ylab("") + xlab("") + ggtitle("by LG thermometer rating") +
  theme(axis.text.y=element_text(color="black", size=9),
        axis.text.x=element_text(color="black", size=10, hjust=c(0,1)),
        plot.title=element_text(color="black", size=13, hjust=.5),
        axis.ticks=element_line(size=.5, color="gray55"))
Figure2g <- Figure2g + coord_cartesian(clip = "off") + draw_label("Effect of two-person relationship", fontface='plain', size=13, angle=90, x=-.23, y=0)
Figure2h <- ggplot(fd_serv_ideo, aes(x=x, y=est, ymin=lo, ymax=hi, group=iv)) +
  geom_hline(yintercept=0, color="gray55", linetype=1, size=.5) +
  geom_ribbon(alpha=.2) +
  geom_line(size=1) +
  scale_y_continuous(limits=c(-.3,.3), breaks=c(-.3, -.2, -.1, 0,.1,.2,.3), minor_breaks=NULL) +
  scale_x_continuous(limits=c(0,1), breaks=c(0,1), labels=c("Very\nliberal","Very\nconservative"), minor_breaks=c(.25,.5,.75)) +
  ylab("") + xlab("") + ggtitle("by ideology") +
  theme(axis.text.y=element_text(color="black", size=9), axis.ticks.y=element_blank(),
        axis.text.x=element_text(color="black", size=10, hjust=c(0,1)),
        plot.title=element_text(color="black", size=13, hjust=.5),
        axis.ticks=element_line(size=.5, color="gray55"))
Figure2i <- ggplot(fd_serv_pid7, aes(x=x, y=est, ymin=lo, ymax=hi, group=iv)) +
  geom_hline(yintercept=0, color="gray55", linetype=1, size=.5) +
  geom_ribbon(alpha=.2) +
  geom_line(size=1) +
  scale_y_continuous(limits=c(-.3,.3), breaks=c(-.3, -.2, -.1, 0,.1,.2,.3), minor_breaks=NULL) +
  scale_x_continuous(limits=c(0,1), breaks=c(0,1), labels=c("Strong\nDemocrat","Strong\nRepublican"), minor_breaks=c(.25,.5,.75)) +
  ylab("") + xlab("") + ggtitle("by party identity") +
  theme(axis.text.y=element_text(color="black", size=9), axis.ticks.y=element_blank(),
        axis.text.x=element_text(color="black", size=10, hjust=c(0,1)),
        plot.title=element_text(color="black", size=13, hjust=.5),
        axis.ticks=element_line(size=.5, color="gray55"))


title1 <- ggdraw() + draw_label("\nEffect on feeling similar to the gay men", fontface='bold', size=17)
title2 <- ggdraw() + draw_label("\nEffect on support for the gay men's case", fontface='bold', size=17)
title3 <- ggdraw() + draw_label("\nEffect on support for required service", fontface='bold', size=17)

p1 <- plot_grid(Figure2a, Figure2b, Figure2c, nrow=1, ncol=3)
p2 <- plot_grid(Figure2d, Figure2e, Figure2f, nrow=1, ncol=3)
p3 <- plot_grid(Figure2g, Figure2h, Figure2i, nrow=1, ncol=3)

plot_grid(title1, p1, title2, p2, title3, p3, ncol=1, rel_heights=c(0.15, 1, 0.15, 1, 0.15, 1)) # rel_heights values control title margins
ggsave("Figure2.pdf", width=10.5, height=20)







##########################################################################################
# STUDY 2: AMCEs
##########################################################################################

# Conjoint analysis
library(cregg)

# AMCEs on feeling similar to teacher
f2 <- cj(Study2.s, similar.teacher~ex_open+ex_poly+ex_gender+ex_race+ex_years+ex_school)
f2$color <- ifelse(f2$estimate==0, "A", "B")
f2$level <- factor(f2$level, levels=rev(levels(f2$level)))

Figure3a <- ggplot(f2, aes(x=level, y=estimate, ymin=lower, ymax=upper, fill=color)) +
  geom_hline(yintercept=0, color="gray55", linetype=1, size=.5) +
  geom_errorbar(position=position_dodge(width=.45), width=.1)+
  geom_point(position=position_dodge(width=.45), size=3.5, shape=21) +
  scale_fill_manual(values=c("white", "gray30")) +
  scale_y_continuous(limits=c(-.3,.3), breaks=c(-.3, -.2, -.1, 0,.1,.2,.3), minor_breaks=NULL) + 
  ggtitle("Effect on feeling similar to the gay teacher\n") +
  ylab("\nAverage Marginal Component Effect") + xlab("") +
  coord_flip() +
  theme(axis.title=element_text(size=11),
        plot.title=element_text(size=12, hjust=.5, face="bold"),
        axis.ticks=element_line(size=.5, color="gray55"),
        axis.text.x=element_text(size=9, color="black"),
        axis.text.y=element_text(size=10, color="black"),
        legend.position="none")
Figure3a <- Figure3a + facet_grid(feature~., scales="free", space="free_y", margins=F) + theme(strip.text=element_blank(), panel.spacing=unit(.4, "lines"))

# AMCEs on supporting teacher in case
f1 <- cj(Study2.s, case.outcome~ex_open+ex_poly+ex_gender+ex_race+ex_years+ex_school)
f1$color <- ifelse(f1$estimate==0, "A", "B")
f1$level <- factor(f1$level, levels=rev(levels(f1$level)))

Figure3b <- ggplot(f1, aes(x=level, y=estimate, ymin=lower, ymax=upper, fill=color)) +
	geom_hline(yintercept=0, color="gray55", linetype=1, size=.5) +
	geom_errorbar(position=position_dodge(width=.45), width=.1)+
	geom_point(position=position_dodge(width=.45), size=3.5, shape=21) +
	scale_fill_manual(values=c("white", "gray30")) +
	scale_y_continuous(limits=c(-.3,.3), breaks=c(-.3, -.2, -.1, 0,.1,.2,.3), minor_breaks=NULL) + 
	ggtitle("Effect on support for the gay teacher's case\n") +
	ylab("\nAverage Marginal Component Effect") + xlab("") +
	coord_flip() +
		theme(axis.title=element_text(size=11),
		plot.title=element_text(size=12, hjust=.5, face="bold"),
		axis.ticks=element_line(size=.5, color="gray55"),
		axis.text.x=element_text(size=9, color="black"),
		axis.text.y=element_text(size=10, color="black"),
		legend.position="none")
Figure3b <- Figure3b + facet_grid(feature~., scales="free", space="free_y", margins=F) + theme(strip.text=element_blank(), panel.spacing=unit(.4, "lines"))


# AMCEs on LGB job discrimination
f3 <- cj(Study2.s, lgb.jobs~ex_open+ex_poly+ex_gender+ex_race+ex_years+ex_school)
f3$color <- ifelse(f3$estimate==0, "A", "B")
f3$level <- factor(f3$level, levels=rev(levels(f3$level)))

Figure3c <- ggplot(f3, aes(x=level, y=estimate, ymin=lower, ymax=upper, fill=color)) +
	geom_hline(yintercept=0, color="gray55", linetype=1, size=.5) +
	geom_errorbar(position=position_dodge(width=.45), width=.1)+
	geom_point(position=position_dodge(width=.45), size=3.5, shape=21) +
	scale_fill_manual(values=c("white", "gray30")) +
	scale_y_continuous(limits=c(-.3,.3), breaks=c(-.3, -.2, -.1, 0,.1,.2,.3), minor_breaks=NULL) + 
	ggtitle("Effect on support for LGB job protections\n") +
	ylab("\nAverage Marginal Component Effect") + xlab("") +
	coord_flip() +
		theme(axis.title=element_text(size=11),
		plot.title=element_text(size=12, hjust=.5, face="bold"),
		axis.ticks=element_line(size=.5, color="gray55"),
		axis.text.x=element_text(size=9, color="black"),
		axis.text.y=element_text(size=10, color="black"),
		legend.position="none")
Figure3c <- Figure3c + facet_grid(feature~., scales="free", space="free_y", margins=F) + theme(strip.text=element_blank(), panel.spacing=unit(.4, "lines"))


(Figure3a + Figure3b + Figure3c)
ggsave("Figure3.pdf", width=18, height=8)







##########################################################################################
# STUDY 2: Do other attributes of LGB people moderate these effects?
##########################################################################################

# AIME on feeling similar to gay teacher
aime.similar <- CausalANOVA(formula=similar.teacher~ex_poly+ex_gender+ex_race+ex_years+ex_school+ex_open,
                  int2.formula=~ex_poly:ex_gender+ex_poly:ex_race+ex_poly:ex_years+ex_poly:ex_school+
                    ex_open:ex_gender+ex_open:ex_race+ex_open:ex_years+ex_open:ex_school,
                  family="gaussian", nway=2, screen=F, collapse=F,
                  data=Study2.s)

# ACE of committed relationship on feeling similar to gay teacher
ce <- ConditionalEffect(aime.similar, treat.fac="ex_open", cond.fac="ex_race")
ce_similar_open <- rbind(ce[[1]][[1]][2,], ce[[1]][[2]][2,], ce[[1]][[3]][2,])
ce <- ConditionalEffect(aime.similar, treat.fac="ex_open", cond.fac="ex_gender")
ce_similar_open <- rbind(ce_similar_open, ce[[1]][[1]][2,], ce[[1]][[2]][2,])
ce <- ConditionalEffect(aime.similar, treat.fac="ex_open", cond.fac="ex_years")
ce_similar_open <- rbind(ce_similar_open, ce[[1]][[1]][2,], ce[[1]][[2]][2,], ce[[1]][[3]][2,])
ce <- ConditionalEffect(aime.similar, treat.fac="ex_open", cond.fac="ex_school")
ce_similar_open <- rbind(ce_similar_open, ce[[1]][[1]][2,], ce[[1]][[2]][2,], ce[[1]][[3]][2,])
ce_similar_open <- as.data.frame(ce_similar_open)
ce_similar_open$group <- c(rep("2ex_race",3), rep("1ex_gender",2), rep("3ex_years",3), rep("4ex_school",3))
names(ce_similar_open) <- c("ce", "sd", "lower", "upper", "group")
ce_similar_open$cond <- c("is White", "is Black", "is Latinx", "is a woman", "is a man", "is in a 3 year\nrelationship", "is in a 6 year\nrelationship", "is in a 10 year\nrelationship", "worked at an\nelementary school", "worked at a\nmiddle school", "worked at a\nhigh school") 
ce_similar_open$cond <- factor(ce_similar_open$cond, levels=c("is Latinx", "is Black", "is White",  "is a man", "is a woman", "is in a 10 year\nrelationship", "is in a 6 year\nrelationship", "is in a 3 year\nrelationship", "worked at a\nhigh school", "worked at a\nmiddle school", "worked at an\nelementary school"))
# Plot ACEs of committed relationship on feeling similar to gay teacher
g_similar_open <- ggplot(ce_similar_open, aes(x=cond, y=ce, ymin=lower, ymax=upper)) +
  geom_hline(yintercept=0, color="gray55", linetype=1, size=.5) +
  geom_errorbar(position=position_dodge(width=.45), width=.1)+
  geom_point(position=position_dodge(width=.45), size=3, shape=21, fill="black") +
  scale_y_continuous(limits=c(-.3,.3), breaks=c(-.3, -.2, -.1, 0,.1,.2,.3), minor_breaks=NULL) + 
  ylab("\nEffect of committed relationship\non feeling similar to the gay teacher") + xlab("") +
  ggtitle("(a)") +
  coord_flip() +
  theme(axis.title=element_text(size=11),
        plot.title=element_text(size=13, hjust=.5),
        axis.ticks=element_line(size=.5, color="gray55"),
        axis.text.x=element_text(size=9, color="black"),
        axis.text.y=element_text(size=10, color="black"))
g_similar_open <- g_similar_open + facet_grid(group~., scales="free", space="free_y", margins=F) + theme(strip.text=element_blank(), panel.spacing=unit(1, "lines"))

# ACE of two-person relationship on feeling similar to gay teacher
ce <- ConditionalEffect(aime.similar, treat.fac="ex_poly", cond.fac="ex_race")
ce_similar_poly <- rbind(ce[[1]][[1]][2,], ce[[1]][[2]][2,], ce[[1]][[3]][2,])
ce <- ConditionalEffect(aime.similar, treat.fac="ex_poly", cond.fac="ex_gender")
ce_similar_poly <- rbind(ce_similar_poly, ce[[1]][[1]][2,], ce[[1]][[2]][2,])
ce <- ConditionalEffect(aime.similar, treat.fac="ex_poly", cond.fac="ex_years")
ce_similar_poly <- rbind(ce_similar_poly, ce[[1]][[1]][2,], ce[[1]][[2]][2,], ce[[1]][[3]][2,])
ce <- ConditionalEffect(aime.similar, treat.fac="ex_poly", cond.fac="ex_school")
ce_similar_poly <- rbind(ce_similar_poly, ce[[1]][[1]][2,], ce[[1]][[2]][2,], ce[[1]][[3]][2,])
ce_similar_poly <- as.data.frame(ce_similar_poly)
ce_similar_poly$group <- c(rep("2ex_race",3), rep("1ex_gender",2), rep("3ex_years",3), rep("4ex_school",3))
names(ce_similar_poly) <- c("ce", "sd", "lower", "upper", "group")
ce_similar_poly$cond <- c("is White", "is Black", "is Latinx", "is a woman", "is a man", "is in a 3 year\nrelationship", "is in a 6 year\nrelationship", "is in a 10 year\nrelationship", "worked at an\nelementary school", "worked at a\nmiddle school", "worked at a\nhigh school") 
ce_similar_poly$cond <- factor(ce_similar_poly$cond, levels=c("is Latinx", "is Black", "is White",  "is a man", "is a woman", "is in a 10 year\nrelationship", "is in a 6 year\nrelationship", "is in a 3 year\nrelationship", "worked at a\nhigh school", "worked at a\nmiddle school", "worked at an\nelementary school"))
# Plot ACEs of two-person relationship on feeling similar to gay teacher
g_similar_poly <- ggplot(ce_similar_poly, aes(x=cond, y=ce, ymin=lower, ymax=upper)) +
  geom_hline(yintercept=0, color="gray55", linetype=1, size=.5) +
  geom_errorbar(position=position_dodge(width=.45), width=.1)+
  geom_point(position=position_dodge(width=.45), size=3, shape=23, fill="black") +
  scale_y_continuous(limits=c(-.3,.3), breaks=c(-.3, -.2, -.1, 0,.1,.2,.3), minor_breaks=NULL) + 
  ylab("\nEffect of two-person relationship\non feeling similar to the gay teacher") + xlab("") +
  ggtitle("(d)") +
  coord_flip() +
  theme(axis.title=element_text(size=11),
        plot.title=element_text(size=13, hjust=.5),
        axis.ticks=element_line(size=.5, color="gray55"),
        axis.text.x=element_text(size=9, color="black"),
        axis.text.y=element_text(size=10, color="black"))
g_similar_poly <- g_similar_poly + facet_grid(group~., scales="free", space="free_y", margins=F) + theme(strip.text=element_blank(), panel.spacing=unit(1, "lines"))



# AIME on support for gay teacher
aime.outcome <- CausalANOVA(formula=case.outcome~ex_poly+ex_gender+ex_race+ex_years+ex_school+ex_open,
                            int2.formula=~ex_poly:ex_gender+ex_poly:ex_race+ex_poly:ex_years+ex_poly:ex_school+
                              ex_open:ex_gender+ex_open:ex_race+ex_open:ex_years+ex_open:ex_school,
                            family="gaussian", nway=2, screen=F, collapse=F,
                            data=Study2.s)

# ACE of committed relationship on support for gay teacher
ce <- ConditionalEffect(aime.outcome, treat.fac="ex_open", cond.fac="ex_race")
ce_outcome_open <- rbind(ce[[1]][[1]][2,], ce[[1]][[2]][2,], ce[[1]][[3]][2,])
ce <- ConditionalEffect(aime.outcome, treat.fac="ex_open", cond.fac="ex_gender")
ce_outcome_open <- rbind(ce_outcome_open, ce[[1]][[1]][2,], ce[[1]][[2]][2,])
ce <- ConditionalEffect(aime.outcome, treat.fac="ex_open", cond.fac="ex_years")
ce_outcome_open <- rbind(ce_outcome_open, ce[[1]][[1]][2,], ce[[1]][[2]][2,], ce[[1]][[3]][2,])
ce <- ConditionalEffect(aime.outcome, treat.fac="ex_open", cond.fac="ex_school")
ce_outcome_open <- rbind(ce_outcome_open, ce[[1]][[1]][2,], ce[[1]][[2]][2,], ce[[1]][[3]][2,])
ce_outcome_open <- as.data.frame(ce_outcome_open)
ce_outcome_open$group <- c(rep("2ex_race",3), rep("1ex_gender",2), rep("3ex_years",3), rep("4ex_school",3))
names(ce_outcome_open) <- c("ce", "sd", "lower", "upper", "group")
ce_outcome_open$cond <- c("is White", "is Black", "is Latinx", "is a woman", "is a man", "is in a 3 year\nrelationship", "is in a 6 year\nrelationship", "is in a 10 year\nrelationship", "worked at an\nelementary school", "worked at a\nmiddle school", "worked at a\nhigh school") 
ce_outcome_open$cond <- factor(ce_outcome_open$cond, levels=c("is Latinx", "is Black", "is White",  "is a man", "is a woman", "is in a 10 year\nrelationship", "is in a 6 year\nrelationship", "is in a 3 year\nrelationship", "worked at a\nhigh school", "worked at a\nmiddle school", "worked at an\nelementary school"))
# Plot ACEs of committed relationship on support for gay teacher
g_outcome_open <- ggplot(ce_outcome_open, aes(x=cond, y=ce, ymin=lower, ymax=upper)) +
  geom_hline(yintercept=0, color="gray55", linetype=1, size=.5) +
  geom_errorbar(position=position_dodge(width=.45), width=.1)+
  geom_point(position=position_dodge(width=.45), size=3, shape=21, fill="black") +
  scale_y_continuous(limits=c(-.3,.3), breaks=c(-.3, -.2, -.1, 0,.1,.2,.3), minor_breaks=NULL) + 
  ylab("\nEffect of committed relationship\non support for the gay teacher's case") + xlab("") +
  ggtitle("(b)") +
  coord_flip() +
  theme(axis.title=element_text(size=11),
        plot.title=element_text(size=13, hjust=.5),
        axis.ticks=element_line(size=.5, color="gray55"),
        axis.text.x=element_text(size=9, color="black"),
        axis.text.y=element_text(size=10, color="black"))
g_outcome_open <- g_outcome_open + facet_grid(group~., scales="free", space="free_y", margins=F) + theme(strip.text=element_blank(), panel.spacing=unit(1, "lines"))


# ACE of two-person relationship on support for gay teacher
ce <- ConditionalEffect(aime.outcome, treat.fac="ex_poly", cond.fac="ex_race")
ce_outcome_poly <- rbind(ce[[1]][[1]][2,], ce[[1]][[2]][2,], ce[[1]][[3]][2,])
ce <- ConditionalEffect(aime.outcome, treat.fac="ex_poly", cond.fac="ex_gender")
ce_outcome_poly <- rbind(ce_outcome_poly, ce[[1]][[1]][2,], ce[[1]][[2]][2,])
ce <- ConditionalEffect(aime.outcome, treat.fac="ex_poly", cond.fac="ex_years")
ce_outcome_poly <- rbind(ce_outcome_poly, ce[[1]][[1]][2,], ce[[1]][[2]][2,], ce[[1]][[3]][2,])
ce <- ConditionalEffect(aime.outcome, treat.fac="ex_poly", cond.fac="ex_school")
ce_outcome_poly <- rbind(ce_outcome_poly, ce[[1]][[1]][2,], ce[[1]][[2]][2,], ce[[1]][[3]][2,])
ce_outcome_poly <- as.data.frame(ce_outcome_poly)
ce_outcome_poly$group <- c(rep("2ex_race",3), rep("1ex_gender",2), rep("3ex_years",3), rep("4ex_school",3))
names(ce_outcome_poly) <- c("ce", "sd", "lower", "upper", "group")
ce_outcome_poly$cond <- c("is White", "is Black", "is Latinx", "is a woman", "is a man", "is in a 3 year\nrelationship", "is in a 6 year\nrelationship", "is in a 10 year\nrelationship", "worked at an\nelementary school", "worked at a\nmiddle school", "worked at a\nhigh school") 
ce_outcome_poly$cond <- factor(ce_outcome_poly$cond, levels=c("is Latinx", "is Black", "is White",  "is a man", "is a woman", "is in a 10 year\nrelationship", "is in a 6 year\nrelationship", "is in a 3 year\nrelationship", "worked at a\nhigh school", "worked at a\nmiddle school", "worked at an\nelementary school"))
# Plot ACEs of two-person relationship on support for gay teacher
g_outcome_poly <- ggplot(ce_outcome_poly, aes(x=cond, y=ce, ymin=lower, ymax=upper)) +
  geom_hline(yintercept=0, color="gray55", linetype=1, size=.5) +
  geom_errorbar(position=position_dodge(width=.45), width=.1)+
  geom_point(position=position_dodge(width=.45), size=3, shape=23, fill="black") +
  scale_y_continuous(limits=c(-.3,.3), breaks=c(-.3, -.2, -.1, 0,.1,.2,.3), minor_breaks=NULL) + 
  ylab("\nEffect of two-person relationship\non support for the gay teacher's case") + xlab("") +
  ggtitle("(e)") +
  coord_flip() +
  theme(axis.title=element_text(size=11),
        plot.title=element_text(size=13, hjust=.5),
        axis.ticks=element_line(size=.5, color="gray55"),
        axis.text.x=element_text(size=9, color="black"),
        axis.text.y=element_text(size=10, color="black"))
g_outcome_poly <- g_outcome_poly + facet_grid(group~., scales="free", space="free_y", margins=F) + theme(strip.text=element_blank(), panel.spacing=unit(1, "lines"))



# AIME on LGB job protections
aime.jobs <- CausalANOVA(formula=lgb.jobs~ex_poly+ex_gender+ex_race+ex_years+ex_school+ex_open,
                            int2.formula=~ex_poly:ex_gender+ex_poly:ex_race+ex_poly:ex_years+ex_poly:ex_school+
                              ex_open:ex_gender+ex_open:ex_race+ex_open:ex_years+ex_open:ex_school,
                            family="gaussian", nway=2, screen=F, collapse=F,
                            data=Study2.s)

# ACE of committed relationship on support for lgb job protections
ce <- ConditionalEffect(aime.jobs, treat.fac="ex_open", cond.fac="ex_race")
ce_jobs_open <- rbind(ce[[1]][[1]][2,], ce[[1]][[2]][2,], ce[[1]][[3]][2,])
ce <- ConditionalEffect(aime.jobs, treat.fac="ex_open", cond.fac="ex_gender")
ce_jobs_open <- rbind(ce_jobs_open, ce[[1]][[1]][2,], ce[[1]][[2]][2,])
ce <- ConditionalEffect(aime.jobs, treat.fac="ex_open", cond.fac="ex_years")
ce_jobs_open <- rbind(ce_jobs_open, ce[[1]][[1]][2,], ce[[1]][[2]][2,], ce[[1]][[3]][2,])
ce <- ConditionalEffect(aime.jobs, treat.fac="ex_open", cond.fac="ex_school")
ce_jobs_open <- rbind(ce_jobs_open, ce[[1]][[1]][2,], ce[[1]][[2]][2,], ce[[1]][[3]][2,])
ce_jobs_open <- as.data.frame(ce_jobs_open)
ce_jobs_open$group <- c(rep("2ex_race",3), rep("1ex_gender",2), rep("3ex_years",3), rep("4ex_school",3))
names(ce_jobs_open) <- c("ce", "sd", "lower", "upper", "group")
ce_jobs_open$cond <- c("is White", "is Black", "is Latinx", "is a woman", "is a man", "is in a 3 year\nrelationship", "is in a 6 year\nrelationship", "is in a 10 year\nrelationship", "worked at an\nelementary school", "worked at a\nmiddle school", "worked at a\nhigh school") 
ce_jobs_open$cond <- factor(ce_jobs_open$cond, levels=c("is Latinx", "is Black", "is White",  "is a man", "is a woman", "is in a 10 year\nrelationship", "is in a 6 year\nrelationship", "is in a 3 year\nrelationship", "worked at a\nhigh school", "worked at a\nmiddle school", "worked at an\nelementary school"))
# Plot ACEs of committed relationship on support for lgb job protections
g_jobs_open <- ggplot(ce_jobs_open, aes(x=cond, y=ce, ymin=lower, ymax=upper)) +
  geom_hline(yintercept=0, color="gray55", linetype=1, size=.5) +
  geom_errorbar(position=position_dodge(width=.45), width=.1)+
  geom_point(position=position_dodge(width=.45), size=3, shape=21, fill="black") +
  scale_y_continuous(limits=c(-.3,.3), breaks=c(-.3, -.2, -.1, 0,.1,.2,.3), minor_breaks=NULL) + 
  ylab("\nEffect of committed relationship\non support for LGB job protections") + xlab("") +
  ggtitle("(c)") +
  coord_flip() +
  theme(axis.title=element_text(size=11),
        plot.title=element_text(size=13, hjust=.5),
        axis.ticks=element_line(size=.5, color="gray55"),
        axis.text.x=element_text(size=9, color="black"),
        axis.text.y=element_text(size=10, color="black"))
g_jobs_open <- g_jobs_open + facet_grid(group~., scales="free", space="free_y", margins=F) + theme(strip.text=element_blank(), panel.spacing=unit(1, "lines"))


# ACE of two-person relationship on support for lgb job protections
ce <- ConditionalEffect(aime.jobs, treat.fac="ex_poly", cond.fac="ex_race")
ce_jobs_poly <- rbind(ce[[1]][[1]][2,], ce[[1]][[2]][2,], ce[[1]][[3]][2,])
ce <- ConditionalEffect(aime.jobs, treat.fac="ex_poly", cond.fac="ex_gender")
ce_jobs_poly <- rbind(ce_jobs_poly, ce[[1]][[1]][2,], ce[[1]][[2]][2,])
ce <- ConditionalEffect(aime.jobs, treat.fac="ex_poly", cond.fac="ex_years")
ce_jobs_poly <- rbind(ce_jobs_poly, ce[[1]][[1]][2,], ce[[1]][[2]][2,], ce[[1]][[3]][2,])
ce <- ConditionalEffect(aime.jobs, treat.fac="ex_poly", cond.fac="ex_school")
ce_jobs_poly <- rbind(ce_jobs_poly, ce[[1]][[1]][2,], ce[[1]][[2]][2,], ce[[1]][[3]][2,])
ce_jobs_poly <- as.data.frame(ce_jobs_poly)
ce_jobs_poly$group <- c(rep("2ex_race",3), rep("1ex_gender",2), rep("3ex_years",3), rep("4ex_school",3))
names(ce_jobs_poly) <- c("ce", "sd", "lower", "upper", "group")
ce_jobs_poly$cond <- c("is White", "is Black", "is Latinx", "is a woman", "is a man", "is in a 3 year\nrelationship", "is in a 6 year\nrelationship", "is in a 10 year\nrelationship", "worked at an\nelementary school", "worked at a\nmiddle school", "worked at a\nhigh school") 
ce_jobs_poly$cond <- factor(ce_jobs_poly$cond, levels=c("is Latinx", "is Black", "is White",  "is a man", "is a woman", "is in a 10 year\nrelationship", "is in a 6 year\nrelationship", "is in a 3 year\nrelationship", "worked at a\nhigh school", "worked at a\nmiddle school", "worked at an\nelementary school"))
# Plot ACEs of two-person relationship on support for lgb job protections
g_jobs_poly <- ggplot(ce_jobs_poly, aes(x=cond, y=ce, ymin=lower, ymax=upper)) +
  geom_hline(yintercept=0, color="gray55", linetype=1, size=.5) +
  geom_errorbar(position=position_dodge(width=.45), width=.1)+
  geom_point(position=position_dodge(width=.45), size=3, shape=23, fill="black") +
  scale_y_continuous(limits=c(-.3,.3), breaks=c(-.3, -.2, -.1, 0,.1,.2,.3), minor_breaks=NULL) + 
  ylab("\nEffect of two-person relationship\non support for LGB job protections") + xlab("") +
  ggtitle("(f)") +
  coord_flip() +
  theme(axis.title=element_text(size=11),
        plot.title=element_text(size=13, hjust=.5),
        axis.ticks=element_line(size=.5, color="gray55"),
        axis.text.x=element_text(size=9, color="black"),
        axis.text.y=element_text(size=10, color="black"))
g_jobs_poly <- g_jobs_poly + facet_grid(group~., scales="free", space="free_y", margins=F) + theme(strip.text=element_blank(), panel.spacing=unit(1, "lines"))


title1 <- ggdraw() + draw_label("Conditional\neffect of\ncommitted\nrelationship,\nwhen\nteacher...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n", fontface='bold', size=14, hjust=0)
title2 <- ggdraw() + draw_label("Conditional\neffect of\ntwo-person\nrelationship,\nwhen\nteacher...\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n", fontface='bold', size=14, hjust=0)

r1 <- plot_grid(title1, g_similar_open, g_outcome_open, g_jobs_open, nrow=1, ncol=4, rel_widths=c(.5,1,1,1))
r2 <- plot_grid(title2, g_similar_poly, g_outcome_poly, g_jobs_poly, nrow=1, ncol=4, rel_widths=c(.5,1,1,1))
r1 / plot_spacer() /r2 + plot_layout(ncol=1, nrow=3, heights=c(1,.06,1))
ggsave("Figure4.pdf", width=18, height=13)
